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Many-body dissipative particle dynamics is constructed to exhibit vapour-liquid coexistence, with 
a sharp interface, and a vapour phase of vanishingly small density. In this form, the model is an 
unusual example of a soft-sphere liquid with a potential energy built out of local-density dependent 
one-particle self energies. The application to fluid mechanics problems involving free surfaces is 
illustrated by simulation of a pendant drop. 



I. INTRODUCTION 

Dissipative particle dynamics (DPD) is familiar as a 
method of simulating complex fluids at a coarse grained 
level 0, 5 for example block copolymer polymer melts 
0, Q , and surfactant solutions 0, |(| . DPD has also been 
used for multiphase fluid problems, such asphase sepa- 
ration kinetics in binary liquid mixtures H, H, 13 > droplet 
deformation and rupture in shear fields [lfj , and droplets 
on surfaces under the influence of shear fields [11]. The 
advantage of DPD for these kind of problems lies in the 
simplicity of the underlying algorithm, and the physical 
way in which singular events such as droplet rupture are 
captured. Such considerations also make the method at- 
tractive for free-surface fluid dynamics problems. Exam- 
ples of these include various kinds of wetting, spreading, 
wicking, and capillary problems. To be used for these 
kind of problems, DPD needs to be extended to allow for 
vapour-liquid equilibrium. In this way a free surface will 
arise naturally as a vapour-liquid interface, and such an 
interface will possess the same physics as a clean vapour- 
liquid interface. 

To achieve vapour- liquid coexistence in DPD, for a sin- 
gle component system, requires a van der Waals loop 
in the equation of state (EOS) (pressure-density curve). 
However this presents a fundamental limitation for stan- 
dard DPD since the soft interaction forces used in the 
method invariably lead to a predominantly quadratic 
EOS [12. One way around this is the 'many-body' DPD 
method invented by Pagonabarraga and Frenkel Hi El, 
and also investigated by Trofimov et al 15]. In many- 
body DPD, the amplitude of the soft repulsions is made 
to depend on the local density. In this way one can 
achieve a much wider range of possibilities for the EOS. 

In the present work, many-body DPD is developed to 
exhibit vapour-liquid coexistence, with a sharp interface, 
and a vapour phase of vanishingly small density. The ap- 
proach taken is fundamentally the same approach as used 
m Refs.GJandlll but with a somewhat different inter- 
pretation of the same mathematics. Therefore a general 
theory for many-body DPD is described first, which I ar- 
gue involves a fundamental re-interpretation of the DPD 
interaction potentials. The specific implementation for 
vapour-liquid equilibrium is described next, and finally 
the application to free surface problems is illustrated by 
simulation of a pendant droplet. Another application to 



vapour-liquid phase separation kinetics was described in 
an earlier note 161. 



II. GENERAL THEORY 



Dissipative particle dyna mics (DPD) is basically 
Molecular Dynamics 0,^3, with two key innovations. 
The first, and perhaps the most profound, is the use of 
soft interactions. This stands in contrast to the common 
use of interaction potentials corresponding to particles 
with hard cores, for example Lennard-Jones interactions 
or modified hard-sphere interactions. The second inno- 
vation is the use of a momentum conserving thermostat. 
This allows one to simulate at a well defined tempera- 
ture yet preserve hydrodynamics, and this can be impor- 
tant for some problems such as phase separation kinetics. 
The thermostat described below is the original (Espahol- 
Warren) thermostat E3> although the Lowe- Anderson 
thermostat is perhaps simpler and more efficient |20|. In 
the present paper, the focus is on the equilibrium prop- 
erties of many-body DPD models, and the nature of the 
thermostat is unimportant. 

The particles in DPD have positions and velocities 
Vj, where i = 1 to N runs over the set of particles, moving 
in a simulation box of volume V. They move according 
to the kinematic condition dvijdt = Vj, and Newton's 
second law dvi j dt = F j /m^ where m t is the mass of the 
ith particle. Here 
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(1) 



is the total force acting on the ith particle, comprising a 
possible external force F^ext and forces Fy due to the in- 
teraction between the ith and jth particles. The interac- 
tion forces are decomposed into conservative, dissipative 
and random contributions, 



(2) 



The individual contributions all vanish for particle sepa- 
rations larger than some cutoff interaction range r c , and 



all obey Newton's third law so that Fy 
The conservative force is 

Fg = Aw c {r l} )e i0 



0. 



(3) 
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where r„- = r,- - r,, r„- = |r y |, and ey = ry/ry. The 
weight function wc( r ) vanishes for r > r c , and for sim- 
plicity is taken to decreases linearly with particle separa- 
tion, thus wq(t) = (1 — r/r c ). This force corresponds to 
a total potential energy which is a sum of pair potentials 

tf({r<}) = E 3 ->i0(r«) (4) 

where —<j>'{r) = Awc(r), and thus 0(r) = (A/2)(l — 
r/r c ) 2 for standard DPD. The dissipative and random 
forces are F° = -7100 (ty, ) (v y - • e y -)ey and = 
(TjKf.(rij)^je i:) -. In these 7 and a are amplitudes, u>d(V) 
and wn(r) are additional weight functions also vanishing 
for r > r c , Vy = Vj — Vj, and = £ji is pairwise con- 
tinuous white noise with (£ij(t)) = and (£y = 
(SikSji + SuSjk)S(t — t'). The dissipative and random 
forces act as the above-mentioned thermostat provided 
that the weight functions and amplitudes are chosen to 
obey a fluctuation-dissipation theorem: a 2 = 2jk B T and 
wd = (wr) 2 , where k B T is the desired temperature in 
units of Boltzmann's constant k B [l9j ■ The same weight 
function is used as for the conservative forces (basically 
for historical reasons): tor = wc and wd = (wc) 2 - 

Usually all the particles are assumed to have the same 
mass, and to fix units of mass and length a convenient 
choice is to set = r e = 1. Often the units of energy 
and hence time are fixed by setting k B T — 1, but for 
equilibrium simulations it can be convenient to keep k B T 
as a free parameter. 

The integration of the equations of motion is a non- 
trivial matter since one has to manage the random forces 
somehow. For an integration algorithm, Groot and War- 
ren investigated a version of the velocity- Verlet scheme 
used in molecular dynamics simulations 0, but it was 
later shown by den Otter and Clarke that this is not a real 
improvement over a simple Euler type integration scheme 
|21j. More extensive studies have been undertaken by 
Vattulainen et al \22\. But many the problems are ob- 
viated if the Lowe-Anderson thermostat is used, which 
is based on distinctly different physical ideas [2(|. All 
the simulations described below though were carried out 
with the simple velocity- Verlet like algorithm described 
by Groot and Warren. 

For a single-component DPD fluid, the equation of 
state (EOS) gives the pressure p as a function of the 
density p = N/V. For the soft potential given above, the 
EOS is now well established to be 0] 

p = pk B T + aAp 2 (5) 

where a = 0.101 ± 0.001 is very close to the mean field 
prediction «mf = ?r/30 = 0.1047 (see also below). The 
first term in the EOS is an ideal gas term, and the sec- 
ond term is the excess pressure, which is almost perfectly 
quadratic in the density (there is a very small correction 
of order p 3 ). Note though that a A is not the second virial 
coefficient 0] > so the above EOS breaks down as p — > 0. 
It seems that a quadratic EOS like this is unavoidable 
for soft potentials ^3 ■ This represents the fundamental 



limitation to basic DPD mentioned in the introduction. 
Moreover one has to take A > otherwise the pressure 
diverges negatively at high densities, so one is are re- 
stricted to a strictly positive compressibility dp/ dp > 0. 
In fact, making A < throws the DPD pair potential 
into a formal class of catastrophic potentials for which 
it can be rig orously proved there is no thermodynamic 
limit |12l |23j . The situation is not as grim as it might 
seem though since considerable progress can be made for 
applications by introducing different species of particles, 
and allowing them to be differentiated by their repulsion 
amplitudes thus A — > Ay in Eqs. © and 

For the one-component fluid, an obvious way to get 
around the problem of a quadratic EOS is to make the 
amplitude A in the force law dependent on density some- 
how. Such a scheme has been examined by several work- 
ers ^3 01 j and proves to be a simple extension to DPD. 
This many-body DPD requires only a modest additional 
computational cost, but throws open the possibility to 
simulate systems with an arbitrarily complicated EOS. 
Care must be taken though with density-dependent in- 
teractions |24|. The approach described here introduces 
a local density into the amplitude in the force law. By 
being explicit about the construction of the local density, 
this is a 'safe' way to introduce a density-dependence into 
the interations |25|. 

In many-body DPD, I write 

Fg- = \[Afa) + Afa)] «te(r«) ey, (6) 

for a one-component fluid (Trofimov et al describe a 
multi-component generalisation). A partial amplitude 
A(p) is introduced, depending on a weighted local den- 
sity, which I define for the ith particle to be 

Pi = T, j7 n w p( r ij)- (7) 

The weight function w p (r) vanishes for r > r c and for 
convenience is normalised so that J d 3 rw p (r) — 1, al- 
though in principle the normalisation could be absorbed 
into the definition of A(pA. The discounted self contribu- 
tion i = j in Eq. would only add a constant w p (0) to 
TJj, amounting to a constant shift of the argument in the 
definition of A(7Jj) (see Trofimov et al for a more exten- 
sive discussion on this point). The weighted local density 
is readily computed by an additional sweep through the 
neighbour list, hence there is only a modest additional 
computational overhead. If A (jo) = A, the method re- 
duces exactly to the standard DPD model. 

In mean field theory, it is easy to show that the modi- 
fied force law should give an EOS 

Puf = pk B T + a M FA(p)p 2 (8) 

where 

a MF = (2tt/3) J^drr 3 w c {r), (9) 

(ie ctMF = 7r/30 for the standard choice of wc (?*))• Thus, 
in principle, an arbitrary dependence on density can be 
recovered. 
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This is not the end of the story though. The existence 
of a potential energy [/({i-j}) such that = dU/dri 
requires the forces to obey a 'Maxwell relation' of the 
type dFi/drj = dFj/dri. This is a non-trivial require- 
ment since the particle positions appear both directly in 
Eq. © and indirectly through the definition of the local 
density in Eq. (JJJ). One can show a necessary and suf- 
ficient condition for it to be true is that the derivative 
w' p is proportional to wq , so the two weight functions w p 
and wc are not independent. One can then prove from 
the normalisation condition on w p {r) that 

- w' p (r) = w c (r)/2a M F (10) 

where qmf is defined in Eq. J^J. 

What, then, is the corresponding potential? I find that 

U({r i }) = j: i u(Pi) (11) 

where u(p^) is a self energy depending on the local den- 
sity, such that 

u'(p) = a MF A(p). (12) 

Comparing Eq. (|llfl with Eq. 10} , it is clear that there 
has been a profound shift in perspective, from a potential 
function expressed in terms of soft pair potentials, to one 
expressed in terms of density-dependent self energies. 

There have recently been many discussions on the 
thermodynamic consistency of density-dependent inter- 
actions in the literature 0, 13, 1^3, EH- However, for 
the present formulation all thermodynamic relations are 
valid because the underlying potential U({ri}) is a well 
defined, density-independent function of the particle po- 
sitions. This is important because it means for instance 
the virial equation for the pressure or stress tensor, con- 
structed out of the forces, can be used without change. 

If u(p) is a polynomial in p of order n, it is easy to show 
that w(/5i) expands to a sum over (n+l)-body density- 
independent potentials. For u(p) = q;mf^P, standard 
DPD is recovered. 

I wish to emphasise, contrary to some hints in the lit- 
erature 0, 0|, that {U({ri})) is the internal energy 
and not the excess free energy (here (. . .) is a thermal 
average). It follows from Eqs. (JSJ) and that the 
mean field EOS is pmf — pksT + p 2 u'{p). A stan- 
dard thermodynamic result for the true pressure is p = 
pk^T + p 2 f^ ' (p) where f^ x \p) is the excess free energy 

per particle. This shows that the interpretation u = f^ 
is a mean field approximation, and as such will be spoilt 
by correlation effects. Note also that correlations mean, 
typically, (p, t ) ^ p, and (u(p, t )) ^ w((^)) ^ u{p). Trofi- 
mov et al give results for the mean local density, and 
suggest ways that one might improve the correspondence 
between (pj and p. Here I take a different approach, and 
regard ~p i as a convenient intermediate quantity which is 
used to construct the forces; as such it is not important 
that its average differs from p. In practice, like Trofimov 
et al, I find that the mean field EOS for many-body DPD 



15 




FIG. 1: Data collapse of pressure against density, (a) stan- 
dard DPD model, for p = 1-10 and A = 0-50. Straight 
line is fit to data given by Eq. © in the text, (b) Many- 
body DPD model, for p = 1-10, A = 0-50, Br\ = 0- 
10 and Td = 0.5-1.0. The ordinate is the function F = 
(p — p — aAp 2 — 2aBr^p 3 ) / Br% Plotted this way, the data 
collapse onto quasi-universal straight lines, where the slope 
depends primarily on r^. This can be expressed as a quasi- 
universal EOS in Eq. (HTl in the text. 

can be considerably less accurate compared to standard 
DPD. Thus the method always requires calibration to 
determine the true thermodynamic properties (the ap- 
proach of Trofimov et al can be used to achieve a specific 
EOS). 

III. A SPECIFIC MODEL 

I now describe in more detail a specific application of 
these ideas to set up a DPD model which exhibits vapour- 
liquid coexistence. Before this though, there is one more 
technical point to discuss. 

To stabilise the vapour-liquid interface, it is not suffi- 
cient just to have a van der Waals loop in the EOS; one 
must also give consideration to the ranges of the interac- 
tions. Thus simple many-body DPD with a single range 
may not have a stable interface as discussed by Pago- 
nabarraga and Frenkel [d|. The trick employed here 
is to take the standard DPD model, make the soft pair 
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I± 
c 



0.50 0.65 0.75 0.85 1.00 



4.0(5) 4.1(1) 3.07(5) 2.08(5) 1.29(5) 



TABLE I: Density correction 'constant' for many-body term 
in measured EOS, as a function of rd- There is little significant 
dependence on other parameters A, B and p (see eg Fig.|Tf bl). 
Figures in brackets are estimates of the error in the final digit. 



potential attractive, and add on a repulsive many-body 
contribution with a different range < r c . Furthermore 
I choose the simplest form of the many-body repulsion, 
namely a self energy per particle which is quadratic in 
the local density. 

In terms of force laws, I take the standard DPD model 
as specified in Eq. © with A < 0, and add a many-body 
force law of the form 



B{Pi + p j )w c {r lJ )e, 



(13) 



where B > 0. This is Eq. ® with A(p) = 2Bp. The 
weight function in this is chosen to be wc(r) = (1 — r/rd) 
for r < rd- This means that w p (r) = 15/(2irr d )(l — r/rd) 2 
(normalised for three dimensions) is used to construct 
the local density, and «mf = 7rr^/30 for this particular 
interaction. 

In terms of potentials, this model can be interpreted 
as follows. Define a generalised weight function of range 
R via 



w p (r; R) = 15/(2tt.R 3 )(1 - r/Rf 



(14) 



Then define two local densities p and p, constructed using 
this weight function with R — r c and R — rd respectively. 
The self energy per particle for this specific model can be 
written as 




FIG. 2: (a) Data collapse of pressure against density, for pl < 
p < 8, A < and \A\ = 20-40, B = 25 or 40, and r d = 
0.75. The ordinate is the function F = (p — p — aAp 2 — 
2aBr^p 3 ) / Br\ Straight line is fit to data given by Eq. JTHJ 
in the text, (b) Pressure as a function of density for the 
two selected parameter sets in Table |H] The lines are the 
predictions of the fitted EOS, Eq. 1)180 in the text. 



u = (n/30)Ap+ (nr d / 30) Bp 2 . (15) 

This is at most quadratic in the local densities, and thus 
the model could be written out explicitly in terms of two- 
and three-body interaction potentials. From this, the 
mean field EOS is 

Pmf = pk B T+(Tr/30)(A + 2Br d p)p 2 . (16) 

Thus, with A < and B > 0, this EOS has the potential 
to contain a van der Waals loop. The actual EOS differs 
from this systematically as I now describe. 

IV. SIMULATION RESULTS 

I now explore by simulation the properties of the above 
model. First I examine the actual EOS, then vapour- 
liquid coexistence and the properties of a stable vapour- 
liquid interface, and finally illustrate the potential ap- 
plication of the method with a simple pendant droplet 
simulation. Typical simulations presented here are in 
simulation boxes of size 10 3 (units of r c ). 



A. Equation of state 

For B — 0, the standard DPD model, the simulations 
recover the accepted equation of state Eq. JSJ with very 
small corrections of ~ p 3 . Results are shown in Fig.^a). 

For A > and B > a large number of simulations 
were performed. After some experimentation, the data 
was found to collapse to the following quasi-universal be- 
haviour, 

p = p + aAp 2 + 2aBr d p 2 (p- c). (17) 

where a takes the same value as for standard DPD (and 
thus this expression contains the correct B = limit), 
and c is an empirical correction to the density that ap- 
pears in the many-body term. This should be compared 
to the mean field prediction in Eq. 1)16(1 . I find that c 
depends predominantly on r d according to Table HJ A 
representative sample of the data is shown in Fig.^b). 

For some parameter sets, the temperature was found 
to show strong deviations from the nominal k^T = 1, 
as a result of instabilities in the integration algorithm. 



3.0 




FIG. 3: Density-temperature phase diagrams for the two pa- 
rameter sets in Table EI Shown horizontal are tielines com- 
puted from vapour-liquid interface profiles, at several tem- 
peratures (keeping A and B fixed). The short vertical line on 
the IcbT = 1 tieline is the point where the pressure vanishes 
according to the fitted EOS, Eq. 1181 in the text. 

Results were only kept if the measured k^T lay within 
10% of the nominal value. These problems occur if the 
repulsion amplitudes are too large, or the densities too 
high, or r<2 too small. The integration algorithm is that 
described in Groot and Warren Q , with a time step At = 
0.05 and A = 1/2. The instabilities can be vanquished 
by making At smaller. 

The measured equation of state is therefore quite close 
to the predicted mean field equation of state. The main 
difference is a correction to the density dependence of 
the many-body term. This is expected since the pair 
correlation function g(r) < 1 where the repulsions are 
strongest, and thus (p^ < p. This effect has also been 
checked in simulations by monitoring the mean value of 
the local density, with results similar to those reported 
by Trofimov et al |15| . 

B. Vapour-liquid coexistence 

For vapour-liquid coexistence I set A < and B > 
so that there is a van der Waals loop in the EOS. Phase 
separation is found in a range of densities pv < p < Ph 
where py and pr, are the vapour and liquid coexistence 
densities. 

In principle, integration of the EOS gives the free en- 
ergy density from which predictions can be made about 
py and pl- Unfortunately, the EOS must deviate from 
the above fitted form for p <C 1, therefore the vapour 
phase is inadequately characterised. For applications, 
one is most interested in pr, > 1 in coexistence with a 
very dilute vapour. If this is true it is much easier to use 
the EOS to predict the point where the pressure vanishes 
as an estimate of the coexisting liquid phase density, thus 
p{ph) = 0. Using this, one expects liquid densities in the 
range pr, ~ 5 for — A ~ B ~ 30. From here on I have set 
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FIG. 4: Interfacial density profiles for the two parameter sets 
in Table ITT1 for several values of the temperature (keeping A 
and B fixed): feT = 1 (solid line), k&T = 2 (chained line), 
k B T = 0.5 (dashed line). 



set 


A 


B 


Td 


PL W O 


dp/ dp 


'5' 
'6' 


-40 

-40 


40 

25 


0.75 
0.75 


5.08(1) 0.78(5) 4.95(3) 
6.05(1) 0.66(3) 7.45(4) 


49(2) 
47(2) 



TABLE II: The two parameter sets used in subsequent simula- 
tions. The sets are distinguished by the different values of the 
liquid densities pl- The coexisting vapour density pv < 1, 
so these parameters are suitable for free surface simulations. 
Also shown are the interface width w, surface tension a, and 
compressibility at p = pl estimated from the EOS. All results 
are at r c = ksT — 1. 

the range of the many-body repulsion to = 0.75 as a 
mid-range value determined above. 

Since the above EOS was measured for A > 0, one has 
to be careful to check that the scaling collapse still holds. 
One cannot easily measure the EOS within the phase sep- 
aration region, since it is hard to maintain a stable uni- 
form density. Therefore the EOS has been characterised 
for p > pl- A similar data collapse is found to the pre- 
vious section, as is shown in Fig. |2Ja) . In this case, the 
EOS can be fitted by 

p = p + aAp 2 + 2aBr\(p 3 - cp 2 + d) (18) 

where a = 0.101(1) as before, c = 4.16(2) and d = 18(1). 
The value of c is similar to the value obtained previously 
(c = 3.07(5), TableHJ). There is an additional offset term 
d which is about 10% of the density correction term cp 2 
in the region of interest. 

Although a wider parameter space was explored, I con- 
centrate here on two parameter sets that were selected for 
more detailed work. These parameter sets are given in 
Table |n] (first three columns). Fig. El^b) shows the pre- 
diction of the EOS, Eq. ljT%)l . compared directly against 
the measured pressures, for these two parameter sets. 

For these two parameter sets, the coexisting vapour 
and liquid densities were determined from the vapour- 
liquid interface profile simulations described in the next 
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8.0 0.5 1.0 1.5 2.0 2.§'° 

kT 

FIG. 5: Interfacial tension a (left hand axis) and interface 
width w (right hand axis), as a function of temperature ksT, 
for the two parameter sets in Tab le|n] (keeping A and B fixed): 
set '5' (solid lines) and set '6' (dashed lines). 

section, and are shown as a function of k-gT in Fig. [3] 
Where k-gT ^ 1 in these simulations, the values of A and 
B are left at the values in Table ITT1 in other words A and 
B are regarded as absolute interaction energies. Also 
shown in Fig. [3] are the appropriate solutions of p = 
using the EOS in Eq. fTSjl. 

It is clear that the difference p^ — pv gets smaller as T 
increases, as one approaches the expected vapour-liquid 
critical point. At k^T < 1, Pv *C 1 indicating the vapour 
phase is virtually devoid of particles. At k B T — 1 the 
solution to p = for the EOS gives a good estimate of 
the density of the fluid phase. 

I have also used the EOS to estimate the compressibil- 
ity dp/ dp at p = pl, and the values are shown in TablelTTI 
Although the precise value is not important, the fact that 
dp/ dp 3> 1 at the coexisting fluid density (where p r* 0) 
shows that the fluid phase is relatively incompressible, 
similar to a real liquid. 

C. Vapour-liquid interface 

Simulations of the vapour-liquid interface were under- 
taken, by taking an equilibrated volume of fluid in a pe- 
riodic box, at a density close to the p = limit, and re- 
moving the particles in one half of the box. The system 
was allowed to evolve until an equilibrium density profile 
was obtained. Interface profiles and surface tension val- 
ues were measured as described for fluid-fluid interfaces 
|10|. For measurement of density profiles, it was neces- 
sary to stop the interface drifting over time. This was 
achieved by inserting a thin slab of 'frozen' particles at 
one end of the box. 

Fig. 0] shows the interface profiles obtained this way 
for the two selected parameter sets in Table [n] These 
are shown at several different values of k B T keeping A 
and B fixed. The limiting densities on either side of 
the interface were used to construct the tielines discussed 
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(a) (b) 

FIG. 6: Pendant droplet problem: (a) isosurface cut through 
3d density field at p — pl/2, showing the drop profile, and 
(b) drop radius r as a function of height z, computed by the 
method described in the text. The solid line is the Young- 
Laplace equation prediction, the circles are the measured pro- 
file. The 'frozen' support particles at the top of the drop give 
a ragged edge to the top of the isosurface in (a), and distort 
the measured profile for z > 23 in (b). Parameters correspond 
to set '6h' in Table ITTT1 



in the previous section (Fig. [3J. As the temperature is 
increased, the interfacial width w gets broader, and the 
surface tension a drops. Results for a and w are shown 
in Fig. |SJ The width w was quantified by calculating 
the maximum slope, and normalising to the coexistence 
densities, thus 



max \dp/dz\ 

The surface tension is determined from the standard me- 
chanical definition of the pressure tensor |17| . Note again 
that there is no problem with the many-body origin of 
the force laws. The actual forces enter the calculation in 
exactly the same way as standard DPD. 

Low temperature favours a sharp interface, but if the 
temperature is too low, oscillations develop in the profile 
on the liquid side of the interface. This can be seen most 
clearly in Fig.0]for k B T = 0.5. The system has crossed a 
Fisher- Widom line in the phase diagram, and a freezing 
transition is almost certainly nearby. The relative am- 
plitude of the oscillations can be measured, and they are 
typically 10% of the bulk density at k B T = 0.5, but < 2% 
for ksT > 0.75, at least for the two sets of parameters 
studied here. 

Thus I conclude that the two parameter sets given in 
Table [n] provide for a sharp vapour-liquid interface, at 
k B T = 1, with virtually no particles in the vapour phase. 
They are thus well suited to model free surfaces. Table ITT1 
also contains the measured interfacial properties. 
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set 


9 PL 


AE 1/H(f3) 


0" shape O" exact 


'5a' 


0.022 5.08 


8.80(5) 0.378 0.57(3) 


4.9(3) 4.95(3) 


'6f 


0.025 6.05 


9.00(5) 0.378 0.57(3) 


7.0(5) 7.45(4) 


'6h' 


0.033 6.05 


8.75(5) 0.404 0.52(2) 


7.9(4) 7.45(4) 


'6g' 


0.030 6.05 


8.70(5) 0.431 0.48(4) 


6.6(5) 7.45(4) 



TABLE III: Pendant drop profile data. Interaction param- 
eters are taken from Table HTl according to the first digit of 
'set'. The penultimate column is the surface tension com- 
puted from the drop shape, and the final column is the 'exact' 
surface tension from Table ITT1 computed by integration of the 
pressure tensor through a planar interface. 



D. Pendant droplet simulation 

As an example application, I look at the classic pen- 
dant drop problem. The procedure is very similar to the 
one adopted for the DPD multiphase fluid model jicj . 

To set up the pendant droplet, a volume of fluid at a 
density close to the equilibrium liquid density was equili- 
brated, then replicated to construct a cylindrical column 
with the axis parallel to the z-direction. A 'support' was 
constructed by 'freezing' particles in a thin slice at the 
top of the column. A gravitational body force g was in- 
cluded by adding a constant force per particle directed 
along the z-direction away from the support. When the 
system reaches equilibrium, the liquid forms a pendant 
droplet suspended from the support particles. In equi- 
librium, the drop profile (radius as a function of height) 
was obtained as described below. The whole simulation 
takes a couple of minutes on a modern workstation. The 
droplet contains typically ~ 3000 particles. 

The profile was determined as follows. A 3d mesh was 
introduced with a resolution typically < 0.5r c (higher 
resolution was employed in the z-direction). The local 
particle density in each mesh volume element was com- 
puted by averaging over a period of time. This gives a 
3d density field. The droplet can then be imaged as an 
isosurface or level cut through this density field, and a 
typical result is shown in Fig. EJa) . 

To determine the drop radius as a function of height, 
the density field was divided (or 'segmented') into occu- 
pied and unoccupied cells according to whether p(r) > 
Pl/2 or not. The number of occupied cells at each height 
z was used to compute the cross-sectional area of the 
droplet at that height, and therefore the drop radius as a 
function of z. A typical drop profile is shown in Fig. EI b). 

This indirect procedure to determine the droplet radius 
eliminates two possible artefacts. Firstly it removes the 
blurring of the base of the drop by the interface profile, 
which would otherwise be ~ 0.7-0.8r c . Secondly it elim- 
inates effects due to the variation of fluid density with 
height which might otherwise introduce a systematic er- 
ror, if the mean particle number density as a function of 
height was computed directly. Such a variation of density 
with height is to be expected, since the fluid responds to 
the varying pressure field through the EOS (ie it is still 



a compressible fluid, even if only weakly so). 

The drop profile was analysed by normalising with 
respect to the maximum diameter AE, and comparing 
with a set of precalculated profiles as the Bond number 
(3 = Phgb 2 1 ' o varies (where b is the radius of curvature 
of the base of the droplet). The profiles are calculated 
from the Young-Laplace equation as described in earlier 
work [ljj. From the best-fit (3 value, the surface tension 
can be computed from a — p^gAE 2 /H where H(/3) is a 
dimensionless function computed numerically. 

Table IIIII shows the quantities computed for several 
drops, for both parameter sets, and for several values of 
g. Although the surface tensions determined this way 
are not very precise, they are all consistent with the ac- 
curate values calculated directly from the interfacial pro- 
files. The drop profiles all match the measured profiles 
quite accurately, see for example Fig.ffib). 



V. DISCUSSION 

The model developed here can be discussed in several 
contexts. Firstly, it is a new simulation method for fluid 
mechanics problems involving liquids with free surfaces. 
For example, the above pendant droplet problem is a test 
of the static force balance and the results show that the 
DPD fluid obeys the Young-Laplace equation in a non- 
trivial geometry. One can conclude that this particular 
version of many-body DPD offers a viable route for solv- 
ing capillary problems such as the distribution of liquids 
in porous materials. It is clearly possible to address dy- 
namic force balance situations too, but these will require 
further testing and parametrisation, particularly for the 
notorious problem of contact line dynamics. 

Secondly, now that vapour-liquid equilibrium is 
achieved for a basic soft sphere model, one can 'dress' 
the liquid up in various ways such as making the liquid 
particles into polymers or model amphiphiles. In this 
way, new methods can be constructed to simulate com- 
plex fluids with an implicit solvent. These developments 
are the subject of ongoing investigations and will be re- 
ported separately. 

In a third context though, the re-interpretation of 
many-body DPD as a fluid whose potential energy is built 
out of local-density dependent one-particle self energies 
is quite novel from the point of view of liquid state the- 
ory. Most previous work has concentrated on fixed pair 
potentials with hard cores, and only minor attention has 
been paid to soft potentials or density-dependent pair po- 
tentials. The present work though goes some way beyond 
these existing ideas. 

It has long been recognised that an arbitrary [/({r^}) 
can be expanded as a sum over density-independent one- 
body, two-body (pair potential), etc, terms. Normally 
the one-body terms, or self energies, are harmless con- 
stants which can be discarded, and most of the phenom- 
ena observed for liquids can be captured by truncating 
the expansion at the pair potential level. If one allows the 
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pieces in such an expansion to acquire a density depen- 
dence though, then the one-body self energy is no longer 
necessarily a constant, and it is no longer necessary to 
go to the pair potential level to see interesting physics. 
Many-body DPD as described here is an example of pre- 
cisely this. 

The phase behaviour of the present model is also po- 
tentially very interesting. By analogy with related soft- 
core systems such as the Gaussian core model frl |2i| , 
and models for polymers of various architectures |3Ctl3l| . 
the particles in the original DPD model are expected to 
freeze into a variety of ordered phases at low tempera- 
tures and intermediate densities, with a re-entrant fluid 



phase at high densities. The version of many-body DPD 
presented in this paper is constructed to have a signifi- 
cant vapour- liquid coexistence region, as shown in Fig.|3 
but the low temperature ordered phases are presumably 
still present, as indicated by the presence of oscillations in 
the liquid side of the vapour-liquid interface in Fig. 01 In 
such a case, the collision between the vapour-liquid tran- 
sition and these ordered phases could prove to generate 
rather unusual phase behaviour, and the low tempera- 
ture properties of many-body DPD models may well be 
worth further examination. 

I thank R. D. Groot for many discussions in the early 
stages of development of the many-body DPD model. 
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